
TGraphErrors** GetRelativeDifference(TGraphErrors* gr1, TGraphErrors* gr2);

void SetOffsetGraphs(TGraphErrors* gr1, TGraphErrors* gr2);

int draw_sum(){

    TFile* g_out_file=new TFile("Fitted_voxel_NNC.root","READ");
    TFile* g_out_file_bulk=new TFile("/dybfs2/users/lijj/nHBackEndAnalysis/AcciAnalysis/Fitted_NNC.root","READ");
    cout<<g_out_file_bulk<<"  "<<g_out_file<<endl;

    TGraphErrors* Constant_voxel=(TGraphErrors*)g_out_file->Get("Constant");
    TGraphErrors* NGd_voxel=(TGraphErrors*)g_out_file->Get("N_Gd");
    TGraphErrors* NLS_voxel=(TGraphErrors*)g_out_file->Get("N_LS");

    TGraphErrors* Constant_bulk=(TGraphErrors*)g_out_file_bulk->Get("Constant");
    TGraphErrors* NGd_bulk=(TGraphErrors*)g_out_file_bulk->Get("N_Gd");
    TGraphErrors* NLS_bulk=(TGraphErrors*)g_out_file_bulk->Get("N_LS");

    Constant_voxel->SetMarkerColor(2);
    Constant_voxel->SetLineColor(2);
    Constant_voxel->SetMarkerStyle(21);
    Constant_voxel->SetMarkerSize(0.8);

    Constant_bulk->SetMarkerColor(4);
    Constant_bulk->SetLineColor(4);
    Constant_bulk->SetMarkerStyle(21);
    Constant_bulk->SetMarkerSize(0.8);

    SetOffsetGraphs(Constant_voxel,Constant_bulk);

    TCanvas* canvas_Constant=new TCanvas("Constant","Constant",800,600);
    TLegend* lg_Constant=new TLegend();
    lg_Constant->SetTextSize(0.05);
    lg_Constant->AddEntry(Constant_voxel,"Sum of voxels","lpe");
    lg_Constant->AddEntry(Constant_bulk,"Whole AD","lpe");
    Constant_voxel->Draw("AP");
    Constant_bulk->Draw("P");
    Constant_voxel->GetXaxis()->SetNdivisions(510);
    Constant_voxel->GetYaxis()->SetRangeUser(-200,200);
    Constant_voxel->GetYaxis()->SetTitle("Constant");
    Constant_voxel->GetXaxis()->SetTitle("AD Index");
    Constant_voxel->GetXaxis()->CenterTitle();
    Constant_voxel->GetYaxis()->CenterTitle();
    lg_Constant->Draw("same");
    
    TGraphErrors** NGd_rel=GetRelativeDifference(NGd_voxel,NGd_bulk);
    TCanvas* canvas_NGd=new TCanvas("NGd","NGd",800,600);
    TLegend* lg_NGd=new TLegend();
    lg_NGd->SetTextSize(0.05);
    lg_NGd->AddEntry(NGd_rel[0],"Sum of voxels","lpe");
    lg_NGd->AddEntry(NGd_rel[1],"Whole AD","lpe");
    NGd_rel[0]->Draw("AP");
    NGd_rel[1]->Draw("P");
    NGd_rel[0]->GetXaxis()->SetNdivisions(510);
    NGd_rel[0]->GetYaxis()->SetRangeUser(-60,60);
    NGd_rel[0]->GetYaxis()->SetTitle("Relative difference [%]");
    NGd_rel[0]->GetXaxis()->SetTitle("AD Index");
    NGd_rel[0]->GetXaxis()->CenterTitle();
    NGd_rel[0]->GetYaxis()->CenterTitle();
    lg_NGd->Draw("same");
    
    TGraphErrors** NLS_rel=GetRelativeDifference(NLS_voxel,NLS_bulk);
    TCanvas* canvas_NLS=new TCanvas("NLS","NLS",800,600);
    TLegend* lg_NLS=new TLegend();
    lg_NLS->SetTextSize(0.05);
    lg_NLS->AddEntry(NLS_rel[0],"Sum of voxels","lpe");
    lg_NLS->AddEntry(NLS_rel[1],"Whole AD","lpe");
    NLS_rel[0]->Draw("AP");
    NLS_rel[1]->Draw("P");
    NLS_rel[0]->GetXaxis()->SetNdivisions(510);
    NLS_rel[0]->GetYaxis()->SetRangeUser(-10,10);
    NLS_rel[0]->GetYaxis()->SetTitle("Relative difference [%]");
    NLS_rel[0]->GetXaxis()->SetTitle("AD Index");
    NLS_rel[0]->GetXaxis()->CenterTitle();
    NLS_rel[0]->GetYaxis()->CenterTitle();
    lg_NLS->Draw("same");
    
    canvas_Constant->SaveAs("Comparison/Constant_rel.pdf");
    canvas_NGd->SaveAs("Comparison/NGd_rel.pdf");
    canvas_NLS->SaveAs("Comparison/NLS_rel.pdf");

    return 0;
}

TGraphErrors** GetRelativeDifference(TGraphErrors* gr1, TGraphErrors* gr2){

    double *X1=gr1->GetX();
    double *Y1=gr1->GetY();
    //double Y1_mean=gr1->GetMean(1);

    double *X2=gr2->GetX();
    double *Y2=gr2->GetY();
    //double Y2_mean=gr2->GetMean(1);

    int N1=gr1->GetN();
    int N2=gr2->GetN();

    if(N1!=N2) return NULL;

    double X1_rel[N1];
    double Y1_rel[N1];
    double Y1_rel_err[N1];

    double X2_rel[N1];
    double Y2_rel[N1];
    double Y2_rel_err[N1];

    for(int i=0;i<N1;i++){

        double Y1_mean=(Y1[i]+Y2[i])/2.;
        double Y2_mean=(Y1[i]+Y2[i])/2.;

        Y1_rel[i]=(Y1[i]-Y1_mean)/Y1_mean*100;
        Y2_rel[i]=(Y2[i]-Y2_mean)/Y2_mean*100;

        X1_rel[i]=X1[i]-0.1;
        X2_rel[i]=X1[i]+0.1;
    
        Y1_rel_err[i]=gr1->GetErrorY(i)/Y1_mean*100;
        Y2_rel_err[i]=gr2->GetErrorY(i)/Y2_mean*100;
    }
    TGraphErrors* grrel_1=new TGraphErrors(N1,X1_rel,Y1_rel,0,Y1_rel_err);
    TGraphErrors* grrel_2=new TGraphErrors(N1,X2_rel,Y2_rel,0,Y2_rel_err);

    grrel_1->SetMarkerColor(2);
    grrel_1->SetLineColor(2);
    grrel_1->SetMarkerStyle(21);
    grrel_1->SetMarkerSize(0.8);

    grrel_2->SetMarkerColor(4);
    grrel_2->SetLineColor(4);
    grrel_2->SetMarkerStyle(21);
    grrel_2->SetMarkerSize(0.8);

    TGraphErrors* *grrel=new TGraphErrors*[2];
    grrel[0]=grrel_1;
    grrel[1]=grrel_2;
    return grrel;
}

void SetOffsetGraphs(TGraphErrors* gr1, TGraphErrors* gr2){
   int N1=gr1->GetN(); 
   for(int i=0;i<N1;i++){
    gr1->SetPointX(i,gr1->GetPointX(i)-0.1);
    gr2->SetPointX(i,gr2->GetPointX(i)+0.1);
   }
}
